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INTRODUCTION 

It is well known that the concept of diffusion is associated with random motion 
of particles in space, usually denoted as Brownian motion, see e.g. [1-3]. Diffusion 
is considered normal when the mean squared displacement of the particle during a 
time interval becomes, for sufficiently long intervals, a linear function of it. When this 
linearity breaks down, degenerating in a power law with exponent a > different from 
one, diffusion is referred to be anomalous: slow if < a < 1 , fast if a > 1 ; see e.g. [4]. 

According to the classical approach started by Langevin and known as the Einstein- 
Ornstein-Uhlenbeck theory of Brownian motion, normal diffusion and Brownian motion 
are associated with Langevin equation. More specifically, the classical Langevin equa- 
tion addresses the dynamics of a Brownian particle through Newton's law by incor- 
porating the effect of the Stokes fluid friction and that of thermal fluctuations in the 
vicinity of the particle into a random force with suitably assigned properties. These 
properties are derived from the requirement that the particle velocity asymptotically 
attains a stationary Maxwellian distribution. Over the period of the diffusing particle, 
the random force arising from molecular collisions undergoes such a rapid fluctuations 
that it is approximated well by a white noise. For large time intervals t, it emerges that 
the mean squared displacement becomes proportional to t with the diffusion coefficient 
being a half of the proportionality constant, in the one dimensional case. 

The present e-print is a reproduction of the contribution published in 1996, so it rep- 
resents our knowledge of that early time. Since 1996 many papers have appeared on 
the topic in view of the rapidly developing theory of fractional diffusion processes. The 
corresponding author (FM) intends to submit an up-dated review on the topic, so he is 
grateful to arXiv readers for any comment and suggestion they may have on this e-print. 
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In Sect. 1. we summarize the salient mathematical aspects of the classical Langevin 
equation, showing the exponential time decay of the velocity correlation function and 
the linear long-time behaviour of the mean squared displacement. We also consider 
the generalized version of the Langevin equation introduced by Kubo [2] to account 
for a general retarded effect of the frictional force and the two fluctuation-dissipation 
theorems by Kubo. 

In Sect. 2, following the approach originally started by Widom [5] and Case [6] and 
resumed in [3-4] , we shall consider the modification of the Langevin equation on the basis 
of hydrodynamics, which takes into account the added mass and the Basset-Boussinesq 
retarding force. We shall improve the analysis of the previous authors, interpreting the 
retarding force in the framework of Fractional Calculus and providing the analytical 
expressions of the autocorrelation functions (both for velocity and random force) and 
of the mean squared displacement. 

We shall conclude noting that, for "not heavy" Brownian particles, there is the 
possibility for anomalous diffusion, with a > 1 , in a long time interval, before the 
normal diffusion is established. 



I. THE CLASSICAL AND GENERALIZED LANGEVIN EQUATIONS 

According to the classical Langevin approach the dynamics in one dimension for a 
Brownian particle is described by 

f = v, 

= (1.2) 

where m is the particle mass, X = X(t), V = V(t) are the particle position and velocity, 
and F is the force acting on the particle from molecules of the fluid surrounding the 
Brownian particle. 

The force F may be divided into two parts. The first part is the frictional force 
and is taken to be proportional to the particle velocity i.e. 

F v = --V, (1.3) 
a 

where 1/a is the friction coefficient for unit mass. One usually introduces the mobility 
coefficient as 

\i\=ajm. (1.4) 
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If the Stokes law is assumed for a spherical particle of radius a , see e.g. [7], we have 

1/ fi = 6tt a pf v , (1-5) 

where pf and z/ are the density and the kinematic viscosity of the fluid, respectively. 
Introducing the characteristic parameters 

To :=a 2 /u, X -=P P /Pf, (1-6) 
where p p is the particle density, we obtain 

1 9 m 2 
2x r 9 

The second part of the force, arising from rapid thermal fluctuations, is regarded 
as random, independent of the motion of the particle. This part is called the random 
force and is hereafter denoted by R(t). 

Then (1.2) is written as a stochastic equation as 

^ = --V(t) + -R(t), (1.8) 
at a m 

and it is referred to as the classical Langevin equation. 

It is assumed that the stochastic processes V (t) and R(t) be stationary. This means 
that the respective autocorrelation functions Cy and Cr, 

C v (t , t) := ( V(t ) V(to + t)) = C v (t), t>0, (1.9) 

and 

C R {t , t) := ( R(t ) R(t + t)) = C R (t) , t > , (1.10) 

do not depend on to • Hereafter we will assume to = . 

As a consequence, because of the Wiener-Khintchine theorem [3], the power spec- 
tra or power spectral densities Iv(w) and Ir(u>) (uj G M) are provided by the Fourier 
transforms of the respective autocorrelation functions. We write 

/+oo . -1 p+oo 

C v (t)e- lujt dt, C v (t) = — I v (u)e +lujt duj, (1.11) 

-oo ^ J — oo 

and 

/+oo . 1 p+oo 

C R {t)e~ lu)t dt, C R (t) = — I R {u)e+ lujt dw. (1.12) 

-oo J — oo 
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We assume that any process f(t) be causal, i.e. vanishing for t < , so that the 
Fourier transform f(ui) of f(t) is related to the Laplace transform by the identity 



POO 

f(w) = f(s)\ a=iu , f(s):= e- st f(t)dt, seC 

Jo 



(1.13) 



It is assumed that the random force has zero mean and is uncorrelated to the 
particle velocity at initial time t = ; in other words, 

(R(t)} = 0, (V(0)R(t)) = 0, t>0. (1.14) 

Furthermore, if the Brownian particle has been kept for a sufficiently long time in the 
fluid at (absolute) temperature T, the equipartition law 

m(V 2 (0)} = kT, (1.15) 

where k is the Boltzmann constant, is assumed for the energy distribution. 

It can be shown (see below) that the previous assumptions lead to the following 
relevant results 

C v (t) = ( V 2 (0) > e-t/* 7 = — e-V* 7 , (1.16) 

m 

and 

C R (t) = ^{V\0))S(t) = ^8(t). (1.17) 
a a 

The result (1.16) shows that the velocity autocorrelation function decays exponen- 
tially in time with the decay constant a, while (1.17) means that the power spectrum 
of R(t) is to be white, i.e. independent on frequency, resulting 

mkT , 
I R (u) = I R = • (1.18) 



The two results can be generalized for the so-called generalized Langevin equation 
introduced by Kubo [2], 



dV 
~dt 



rt 1 

- ^(t-T)V{r)dr+—R{t), (1.19) 
Jo m 



where the function ^(t) represents a retarded effect of the frictional force. For this case 
Kubo introduced two fluctuation-dissipation theorems that, using the Laplace trans- 
forms, read respectively 

s + 7(sJ 
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and 

C R (s)=m 2 (V 2 (0))-f(s). (1.21) 

From the comparison between (1.8) and (1.19), we recognize that the classical case 
can be obtained from the generalized one interpreting the convolution in (1.19) in the 
generalized sense (see e.g. [8]) and putting 

7(0 = -*(*) <=► 7(«) = (1-22) 

O" O" 

where 5(t) denotes the delta Dirac distribution. 

In Appendix A we prove the statements (1.20) and (1.21), which thus reduce to 
the classical results (1.16) and (1.17) accounting for (1.22). 

It can be readily shown that the mean squared displacement of a particle, starting 
at the origin at to = , is given by 

(X 2 (t)) = 2 f(t-T)C v (T)dT = 2 f dr x I CV(r)dr. (1.23) 
Jo Jo Jo 

For this it is sufficient to recall that X(t) = f* V(t') dt' , and to use the definition (1.9) 
ofCV(t). 

For the classical case, Cy(£) is provided by (1.16), so that we obtain from (1.23) 

( X 2 (t) ) = 2 ( V 2 (0) ) a [t - a (l - e - */* 7 ) ] . (1.24) 

Introducing the diffusion coefficient D as 

D = a{V 2 (0)} =/ikT, (1.25) 

where we have used (1.4) and (1.15), we can deduce for large times the well-known 
property 

(X 2 (t)) = 2Dt|l + (t/o-y 1 } , as t^oo. (1.26) 

The relationship stated in (1.25), which is called the Einstein relation, provides us with 
a very good basis of experimental verification that Brownian motion is in fact related 
to the thermal motion of molecules. We also note the following relevant results for the 

diffusion coefficient 

D= lim ( X y^ = [ C v (t)dt = C v (0). (1.27) 

t^oo 2 1 Jo 
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2. THE FRACTIONAL LANGEVIN EQUATION 



On the basis of hydrodynamics the equation of motion (1.8) is not all correct since 
it ignores the effects of the added mass and of the retarded viscous force, which are due 
to the acceleration of the particle, see e.g. [5-7]. 

The added mass effect introduces a modification in the L.H.S. of (1.2) in that it 
requires to substitute the mass of the particle with the so-called effective mass, namely 

m -> m e := m+^rrif = m ^1 + ^-^ , (2.1) 

where x = Pp/Pf according to (1.6). As a consequence, in order to do not change the 
mobility coefficient in the Stokes drag, we have to introduce a e such that 

fx := a/m = a e /m e , (2.2) 

namely, recalling (1.6-7), 
where To = a 2 jv. 

With respect to the classical analysis, it turns out that the added mass effect, if 
it were present alone, would be only to lengthen the time scale (a — > a e > a ), slowing 
down the exponential decay for the velocity correlation function (1.16) and for the mean 
square displacement (1.24), but without modifying the value of the diffusion coefficient. 

The retarded viscous force effect is due to an additional term to the Stokes drag, 
which is related to the history of the particle acceleration. This additional drag force, 
proposed independently by Boussinesq [9] and Basset [10] in earlier times, is nowadays 
referred to as the Basset force. As a consequence, the frictional force (1.3-5) is to be 
substituted as follows 

F. = -**ap, V \v(t) +'f d -^0i- dr) . (2.4) 
Using (1.6-7) and requiring the causality of the processes, we can re- write (2.4) as 

(2.5) 



9 

F v = ~—m 
2x 



-V(t) + ^=B(t) 



where 



1 [> dV(T)/dT 1 [> d,V(T)/dT 
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The lower limit of the integral has been written as 0~ to account for the possible 
discontinuity in the velocity particle at t = 0. Basing on the Fractional Calculus recalled 
in Appendix B, we can write, see (B. 11-12) and (B.6), 



B(t)=Dl /2 V(t) = <S>_ 1/2 (t) * V(t), 
where D\J 2 denotes the fractional derivative of order 1/2 and 

£-3/2 £-3/2 



(2.7) 



$-i/ 2 (*) := 



T(-l/2) 2v^F ' 



(2.8) 



Then, adding the random force R(t), the complete Langevin equation (1.2) turns 
out to be 

■— [l + V^Dl /2 ] V(t) + — R(t) . (2.9) 



dV_ 
~~dt 



We agree to refer to (2.9) as to the fractional Langevin equation. 

We recognize that our fractional Langevin equation is a particular case of the 
generalized Langevin equation (1.19) with 



lit) 



1 

CT fi 



s(t) - 



£-3/2 



1 



7(s) = — l + v 7 ^^ 72 



(2.10) 



Consequently, we can use (2.10) to compute the correlations functions Cv(t) , Cn(t) 
starting from their Laplace transforms (1.20-21), respectively. Then, the mean squared 
displacement can be derived from Cyif) according to (1.23). 

Let us first consider the random force. The inversion of the Laplace transform 
Cr(s) yields 

C R (t)=m 2 e (V 2 (0)) 7 (t), (2.11) 

where 7(f) is provided by (2.10). We thus recognize that for our fractional Langevin 
equation the random force cannot be longer represented uniquely by a white noise; an 
additional "fractional" noise is present due to the term t -3 / 2 which, as formerly noted 
by Case [6], is to be interpreted in the generalized sense of tempered distributions [8]. 

Let us now consider the velocity correlation. Using (1.20) and (2.10) it turns out 

<U 2 (0)> <U 2 (0)> 



CV{S) s+ [l + V^aVa]/^ s + y/pfesW + l/ae 
where, because of (2.3) and (1.6), 

r __ 9 9/9/ 



(2.12) 



2 X + 1 2 Pp + p f 
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(2.13) 



We note from (2.13) that < (3 < 9, the limiting cases occurring for x = oo and x = 0, 
respectively. We also recognize that the effect of the Basset force is expected to be 
negligible for j3 — > , i. e. for particles which are sufficiently heavy with respect to the 
fluid (p p > p f ). 

As far as we know, at least in this context, an explicit inversion of the Laplace 
transform (2.12) in terms of elementary functions has not yet been carried out. Widom 
[5] and Case [6] have only provided integral representations of the velocity correlation 
function, from which they have derived the long-time asymptotic behaviour (ex t -3 / 2 ). 
In our notation, applying the asymptotic theorem for s — > to (2.12), we get as t — > oo 



C v (t)^(V 2 (0) 



(±y 3/2 = (v 2 (q)) (±y 3/2 



(2.14) 



The presence of such a long-time tail, pointed out also in [3-4], was first observed by 
Alder and Wainwright [11] in a computer simulation of velocity correlation functions. 

The explicit inversion of (2.12) is hereafter carried out, basing on our previous anal- 
ysis of the original and generalized Basset problems, in the framework of the Fractional 
Calculus and Mittag- Leffler functions [12-14]. For this aim let us recall the following 
Laplace transform pairs 



1 



(S 1 /2_ a+ )( s l/2_ a _) • a+ 



a + E 1/2 (a + Vt) - a_ E 1 / 2 (a_ \fi) 



(2.15) 



(s i/ 2 1 _ a )2 + E V2{aVt) [l + 2a 2 t] +2a v / tA, 



where 



denotes the Mittag- Leffler function of order 1/2. In fact, re- writing (2.12) as 



(2.16) 



00 n n f n/2 2 

E 1/2 (aVt) :=E nn/2 + H = e ° erfc (-^) (2.17) 
n=0 \ / ' ) 



V 2 (0)) 



CV{S) ( S 1 /2_ a+ )( s l/2_ a _) 



(2.18) 



where 



■\TB± f/3-4) 1 / 2 1 
a±= VP^IP_^J if p^ A ^ a± = a = L= if /9 = 4 7 (2.19) 



we easily obtain the required Cy{t) in terms of Mittag- Leffler functions, as pointed out 
in (2.15-16). 
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Furthermore, it can be proved that Cy(t) results for t > a decreasing function, 
completely monotonic, i.e. (—l) n Cy\t) > 0, with the asymptotic behaviour given by 
(2.14), for any physical value of (3 . 

In order to compute (X 2 (t) ), according to (1.23) we have to consider the 2-fold 
primitives of the functions in the R.H.S. of (2.15-16), vanishing at t = 0. In particular, 
the repeated integral for the Mittag-Leffler function turns out 



°^ a n t n/2+2 

/ o^(^) = E FR 2T3) 



n=0 

1 
'a 1 



f l/2 4 t 3/2 

E 1/2 (aVt) - 1 - 2a -= - a 2 1 - -a 3 -= 

V7T 3 V7T 



(2.20) 



The asymptotic behaviour of ( X 2 (t) ) as t — > oo can be easier obtained from its Laplace 
transform for s — > , and reads 



where 



{X 2 (t)} = 2Dt|l + (t/(T e )" 1/2 ]| , as t^oo, 



-D = CV(0) = a e ( F 2 (0) ) = fj,kT . 



(2.21) 



(2.22) 



Note that in the RHS of (2.22) we have used the energy equipartition law (1.15) with 
the effective mass and (2.2). 

The explicit expressions of the velocity autocorrelation function and of the dis- 
placement variance are given in [14-15]. 



CONCLUSIONS 



In this paper we have revisited the Brownian motion on the basis of the fractional 
Langevin equation (2.9), which turns out to be a particular case of the generalized 
Langevin equation (1.19) introduced by Kubo on 1966. 

The importance of our approach is to model the Brownian motion more realistically 
than the usual one based on the classical Langevin equation (1.8), in that it takes into 
account also the retarding effects due to hydrodynamic backflow, i.e. the added mass 
and the Basset memory drag, as pointed out in (2.1) and (2.4), respectively. 

On the basis of the two fluctuation- dissipation theorems (recalled in the Appendix 
A) and of the techniques of the Fractional Calculus (recalled in the Appendix B), we 
have provided the analytical expressions of the correlation functions (both for the ran- 
dom force and the particle velocity) and of the mean squared particle displacement. 
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Consequently, the well-known results of the classical theory of the Brownian motion 
have been properly generalized. 

The random force has been shown to be represented by a superposition of the usual 
white noise with a "fractional" noise, as pointed out in (2.10-11), 

The velocity correlation function Cy{t) exhibits a different behaviour from the 
classical case: it is no longer expressed by a simple exponential but by a combination 
of Mittag-Leffler functions of order 1/2, according to (2.15-19). As a consequence, one 
can derive for Cy(t) a slower decay, proportional to t~ 3 / 2 as t — > oo , which indeed is 
more realistic than the usual exponential one, also in view of numerical simulations. 

Finally, the mean squared displacement has been shown to maintain, for sufficiently 
long times, the linear behaviour which is typical of normal diffusion, with the same 
diffusion coefficient of the classical case, as seen in (2.21-22), i.e. (X 2 (t)) ~2Dt. 
However, the Basset memory force, which is responsible of the algebraic decay of the 
velocity correlation function, induces a retarding effect (oc t 1 / 2 ) in the establishing 
of the linear behaviour, which is relevant when the parameter (3 introduced in (2.13) 
is big enough. From numerical computations this effect is seen to be evident when 

< p p < 2 pf, i.e. for "not heavy" Brownian particles; in these cases one can get a best 
fit in a long time interval with the law (X 2 (t) ) ~ 2 D* t a , with < D* a 0-1 < D and 

1 < a < 2 , which appears as a manifestation of fast anomalous diffusion [14-15]. 
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APPENDIX A 



Let us consider the generalized Langevin equation (1.19), that we write as 



R(t) = m V(t) + 7 (£) * V(t) 



(A.l) 



where • denotes time differentiation and * time convolution. The assumption of station- 
arity for the stochastic processes along with the following hypothesis 



(R(t)) = 0, ( V (0) R(t) ) = , £>0, 
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(A.2) 



allows us to derive, by using the Laplace transforms, the two fluctuation- dissipation 
theorems 

(V 2 (0)) 



C v (s) :=(V(0)V(t) 



s + 7(s) 



and 



C R (s) 



\R{0)R{t)) = m'{V'(p))i{8) 



(A3) 



(A4) 



Our derivation is alternative to the original one by Kubo who used Fourier transforms 
[2] ; furthermore, it appears useful for the treatment of our fractional Langevin equation. 



Multiplying both sides of (A.l) by V(0) and averaging, we obtain 
(V(0)V(t)) +1 (t)*{V(0)V(t)) = 0. 
The application of the Laplace transform to both sides of (A. 5) yields 

s < W) W) )-(V 2 (0))+ 7( S ) ( V(0) V(t) ) = , 

from which we just obtain (A. 3). 

Multiplying both sides of (A.l) by R(0) and averaging, we obtain 

C R (t) := ( R(0) R(t) ) = m 2 [( V(0) V(t) ) + 7 (£) * ( V(0) V(t) ) 

Noting that, by the stationary condition, 

(V(0)V(0)) = 0, (V(0)V(t)) = -(V(0)V(t)), 



the application of the Laplace transform to both sides of (A. 7) yields 

C R (s)=m 2 {s(V(0)V(t))-^(s) [s(V(0)V(t))-(V 2 (0))]} 

Since 



( V(0) V(t) ) = -( V(0) Vit) ) = -s( V(0) V{t) ) + ( V\Q) ) , 



we get 



(A5) 



(A6) 



(A7) 

(A8) 

(A9) 
(A10) 



C R (s) = m 2 {s [-sC v (s) + (V 2 (0))--f(s)C v (s)]+-?(s)(V 2 (0))} , (All) 
from which, accounting for (A. 3), we just obtain (A. 4). 
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APPENDIX B 



Here we recall the essentials of Riemann-Liouville Fractional Calculus basing on 
[16-20], and we interpret the Basset force in terms of a fractional derivative of order 
1/2. 

Usually, the starting point to introduce the Riemann-Liouville fractional calculus 
is the well-known Cauchy's iterated formula, which provides the n-fold primitive of a 
given function /(£) in terms of a single integral. If t > c G R, it reads 



J ...J J f(r)drdT 1 ...dr n . J 



ri)i/(*- T ) n " 1 /( T ) dT ' n = 1 > 2 > 



(n-l) 

Here c denotes the point where the primitive is required to vanish along with its first 
n—1 derivatives. The passage from n G IN to a G 1R + is now quite natural taking into 
account that (n — 1)! = T(n) . Consequently we define 

- the fractional integral of f(t) of order a (with starting point c) 

i:f{t):=-?— /(t-rr-V(r)rfr, a G R + . (B.2) 

For a = we define 1° f(t) = f(t) so that 1° = I where I is the identity operator. The 
choice with c = — oo is originally due to Liouville (1832), while with c = to Riemann 
(1847). 

In order to introduce the notion of fractional derivative of order a , we need to 
consider the possibility to change a — > —a in the r.h.s. of (B.2). While the extension 
in (B.l) from n to a > is quite legitimate, the actual proposal requires some care due 
to the convergence of the integral. 

If a denotes any positive real number in the range n — 1 < a < n with n G IN and 
fit) is a sufficiently well-behaved function, one usually defines 

- the fractional derivative of f(f) of order a (with starting point c) 

d n 1 d n r l f(r) 

c J w dt n c J w r(n - a) dt n J c (t - T ) a+1 - n v ' 

A different formula for the fractional derivative, alternative to (B.3), originally 
introduced by Caputo [14-15], is 

D<*f(t):=I™- a ^—f(t) = —^ / . J . K > dr. (BA) 

c J W c dt n J W r(n - a) J c (t - t )«+!-» v ; 
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We note in general that 



^ ir a m ± ir a ^ m , (5.5) 



unless the function f(t) along with its first n — 1 derivatives vanishes at t = c + . 

For causal functions (i.e. vanishing for t < 0) the choice c = is in order. In this 
case it is convenient to introduce the so-called Gel'fand-Shilov distribution [8] 

*A(t):=^6(f), AGC, (5.6) 

where 6(t) is the unit step Heaviside function and T(A) is the Gamma function. For 
A = — n (n = , 1 , . . .) , $>\(t) reduces to the n-derivative (in the generalized sense) of 
the Dirac delta distribution, 

-i—n—l 

^-n(t):=Y^e(t) = S^(t), n=0,l,... (B.7) 

Assuming that the passage of the n-derivative in (B.3) under integral is legitimate, one 
recognizes that, for n — 1 < a < n , 

n-l 

D% fit) = D% f{t) + f ik) (0 + ) *(*-«+!)(*) , (B.S) 

fc=0 

and, using the (generalized) technique of Laplace transforms, 

n-l 

C {Do fit)} = a a f(a) - ^~ 1 ~ k f (k) (0 + ) ■ (B.9) 

fc=0 

Let us now consider the causal restriction of the Basset force. We easily recognize 
in (2.6) that 

1 f* dV(r)/dr , 1 f* dV(T)/dr , ^ « ,„ , 

Consequently, we can write the following equivalent expressions in terms of the deriva- 
tives of order 1/2 , 

B(t) = bl /2 V(t) + V(0) * 1/2 (t) = D\ /2 V(t) . (B.ll) 
Applying the property (B.9) to (B.ll), we can also write 

B(s) = s l ' 2 V(s) B(t) = $- 1/2 (t) * V(t) , (S.12) 
where the convolution is to be intended in the generalized sense of Gel'fand-Shilov [8]. 
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